SUBROUTINE G_MOT_PHI
!-----------------------------------------------------------------------
!
! Definition von G:  Bewegungs- bzw. Impulsgleichung
!

      use primvar,  only : G, X, XZ, XA, MUp, MD, MR, Gnorm
      use physco,   only : grav, pi, offset, z1, z12
      use config,   only : np, BCflag
      use global,   only : tst
      use geomvar,  only : S_flux, S_vol, S_volA, S_volZ
      use viscvar,  only : Uterm_p, muQ_p
      use advecvar, only : momp_adv
      use RBvar,    only : M_star


      implicit none

      integer :: i
      real(kind=8) :: scale_uphi
      logical, save :: init = .true.

      if (init) then
         write(66,"(a)") "G_mot_phi.f90:   Innenrand: u_phi = Kelper'sch ;  Aussenrand: u_phi = Kepler'sch"
         init = .false.
      end if


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

   ! innere Pseudozellen: i=1 & i=2
      G(MUp,1)        =  X(MUp,1) - sqrt(grav*M_star/X(MR,1))
      Gnorm(MUp,1)    =  max( abs(X(MUp,1)), sqrt(grav*M_star/X(MR,1)) )
      G(MUp,2)        =  X(MUp,2) - sqrt(grav*M_star/X(MR,2))
      Gnorm(MUp,2)    =  max( abs(X(MUp,2)), sqrt(grav*M_star/X(MR,2)) )

   ! aeussere Pseudozellen: i=np, i=np-1 & i=np-2

	   scale_uphi = 1.e-06

       G(MUp,np)       =  X(MUp,np) - sqrt(grav*M_star/X(MR,np)) * scale_uphi
       Gnorm(MUp,np)   =  max( abs(X(MUp,np)), sqrt(grav*M_star/X(MR,np)) * scale_uphi )
       G(MUp,np-1)     =  X(MUp,np-1) - sqrt(grav*M_star/X(MR,np-1)) * scale_uphi 
       Gnorm(MUp,np-1) =  max( abs(X(MUp,np-1)), sqrt(grav*M_star/X(MR,np-1)) * scale_uphi )
       G(MUp,np-2)     =  X(MUp,np-2) - sqrt(grav*M_star/X(MR,np-2)) * scale_uphi
       Gnorm(MUp,np-2) =  max( abs(X(MUp,np-2)), sqrt(grav*M_star/X(MR,np-2)) * scale_uphi )

 !     G(MUp,np)       =  X(MUp,np)
 !     Gnorm(MUp,np)   =  max( abs(X(MUp,np)), sqrt(grav*M_star/X(MR,np)) )
 !     G(MUp,np-1)     =  X(MUp,np-1)
 !     Gnorm(MUp,np-1) =  max( abs(X(MUp,np-1)), sqrt(grav*M_star/X(MR,np-1)) )
 !     G(MUp,np-2)     =  X(MUp,np-2)
!      Gnorm(MUp,np-2) =  max( abs(X(MUp,np-2)), sqrt(grav*M_star/X(MR,np-2)) )

!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

      do i=3,np-3

         G(MUp,i) = S_vol(i)  * X(MUp,i)  * X(MD,i)  * z12*( X(MR,i)  + X(MR,i+1)  )                                 & ! temporal diff.
                  - S_volA(i) * XA(MUp,i) * XA(MD,i) * z12*( XA(MR,i) + XA(MR,i+1) )                                 & ! ...
                  + S_flux(i+1) * XZ(MR,i+1) * momp_adv(i+1) - S_flux(i) * XZ(MR,i) * momp_adv(i)                    & ! advection
                  - z12*pi * (   muQ_p(i+1) * XZ(MR,i+1)**2 * Uterm_p(i+1)   -                                       & ! viscosity
                                 muQ_p(i)   * XZ(MR,i)**2   * Uterm_p(i)     ) * tst                                   ! ...


         Gnorm(MUp,i) = max( offset, &
               abs( S_vol(i)  * X(MUp,i)  * X(MD,i)  * z12*( X(MR,i)  + X(MR,i+1)  ) ),                              & ! temporal diff.
               abs( S_flux(i+1) * XZ(MR,i+1) * momp_adv(i+1) ), abs ( S_flux(i) * XZ(MR,i) * momp_adv(i) ),          & ! advection
               abs( z12*pi * (   muQ_p(i+1) * XZ(MR,i+1)**2 * Uterm_p(i+1)   -                                       & ! viscosity
                                 muQ_p(i)   * XZ(MR,i)**2   * Uterm_p(i)     ) * tst )                               & ! ...
                           )

      end do


END SUBROUTINE G_MOT_PHI

